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Column  vector 

Row  vector 

Kronecker  delta 
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INTRODUCTION 

This  report  is  a  continuation  of  a  paper  by  Kirmser  and 
Wakabayashi  [l  3   on  the  determination  of  satellite  position  vectors. 

In  satellite  photography,    a  number  of  cameras  located  at  dif- 
ferent positions  and  fixed  to  the  earth  yield  plates  which  show  section- 
ally  continuous  satellite  traces.     This  experimental  data  may  then  be 
reduced  to  unit  sighting  vectors. 

The  plates  may  have  simultaneous  tinning  marks,    but  generally 
they  will  not,    although  the  plates  do  overlap  in  tinne.     In  the  latter  case 
it  is  necessary  to  use  sonme  point-matching  technique  to  deternnine 
corresponding  unit  sighting  vectors. 

The  purpose  of  this  report  is  to  determine  the  position  vectors 
of  a  portion  of  the  orbit  of  satellite  i960  Iota  one.    Revolution  No.    121, 
from  experimental  data  that  has  already  been  reduced  to  unit  sighting 
vectors.     The  accuracy  of  the  satellite  locations  and  the  effect  of  small 
variations  in  the  orbit  paranmeters  will  be  examined. 

DETERMINATION  OF  SATELLITE   POSITIONS 

In  general,    simultaneous  sighting  data  from  N  stations  leads  to 
inconsistent  sets  of  equations  (Fig.    1  ).       In  [1  j  it  is  shown  that  the 
inconsistent  sets  of  equations  can  be  solved  in  the  least- square  sense 
by  finding  the  point  S,    such  that  E,    which  is  defined  by  ( l) ,    is  a  minimum. 


Fig.    1.     The  geometl-y  of  simultaneous 
sightings  from  three  stations. 
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E    is  minimized  with  respect  to  S  by  setting -g^ 
This  gives 

{S}   =  A-1     ^   (M.)    {P.   } 


i=l 


where 


A-   I     (M.), 
i=l 


This  theory  has  been  extended  f  1  J    to  include  the  case 
when  the  sighting  data  from  the  N  stations  is  continuous,    or 
consists  of  many  discrete  unsynchronized  observations  (Fig.    Z) . 
Formulas  for    E    and    S    will  now  be  functions  of  Sonne  parameters 
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^  .  ,   j  =   1,    2,    .  .  .^  N  and  the  point- matching  problem  becomes  finding 
a  set  I  .  ,   j  =   1,    2,    .  .  .    ,   N-1  which  minimizes 
N 

E   =  ^  Lp."^  -s'^J  (M.)  {p.  -  S}  =  E(Ii,  ^2.  ....  i^)- 
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Fig.    2.      Point-nnatching  for  continuous 
data  from  N  stations. 


This  leads  to  a  set  of  non-linear  algebraic  equations  in  the 
parameters  4  J  ,    ...    ,    ^^t_,    of  the  form  ' 

j_Pj^  -  s'^J  (M^)   {    Pj^  -  S  }   =  0.       K=   1,  2,    ...    .    N-1  (4) 


Equation  (4)  may  be  solved  using  a  multi- dimensional  Newton's 
method  --a  successive  approximation  technique  in  which  the  initial 
values  of  the  parameters  are  estinnated  and  the  corrections  to  these 
initial  values  are  determined  from  the  matrix  equation 
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In  these  formulas, 
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Because  of  the  iterative  nature  of  this  technique,    it  is  ideally 
suited  for  solving  by  a  computer. 

A  program  was  written  to  solve  the  system  of  non-linear 
algebraic  eq\iations.     In  addition,    the  program,   written  in 
Fortran  II  for  the  IBM  l620,    computed    E    and    S    using  the  adjusted 
parameters.     The  program,    Appendix  I,   was  designed  to  continue 
adjusting  the  parameters  for  a  given  set  of  data  iintil  all  of  the  para- 
meter corrections  became  smaller  than  an  arbitrarily  selected  value. 
This  arbitrary  value  could  be  adjusted  to  improve  the  convergence 
any  time  the  program  was  in  the  computer. 

The  physical  identity  of  the  parameters  is  entirely  arbitrary. 
The  parameter  used  may  be  the  observation  ntimber,    one  co-ordinate 
of  a  sighting  vector,    or  the  time  of  the  observation.     In  fact,    the  para- 
meters used  may  vary  from  station  to  station. 

The  time  of  the  observation  was  chosen  as  the  parameter  with 
the  express  purpose  of  determining  the  timing  error  between  the 
various  stations. 

The  program  was  applied  to  a  set  of  sighting  vectors  from  four 
stations  which  were  located  at  Manhattan,    Kansas;  Oakley,    Kansas; 
St.    Louis,    Missouri;  and  Menlo  Park,    California.     (The  experimental 
data  may  be  foiind  in  Appendix  II.)     The  Manhattan  station  was  the  ref- 
erence station  during  the  adjustment  of  the  parameters.     That  is,    the 
initial  value  of  the  parameter  of  the  Manhattan  station  was  not  adjusted, 
whereas  the  others  were  corrected  to  satisfy  (4).       The  value  of  the 


Manhattan  parameter  was  used  as  the  initial  estimate  of  the  para- 
meters at  the  other  stations. 

As  stated  earlier,  the  program  computed  E  and  S  for  each 
trial.  Since  E  is  the  svtm  of  the  errors  squared,  the  radius  of  an 
error  sphere  nnay  be  determined  by 

Hence  around  each  point  in  space  determined  by    S    there 
is  an  error  sphere  of  radius    dP.       Combining  these  error  spheres 
over  the  time  interval  of  observation  gives  an  irregular  tube  which 
approximates  the  orbit  in  the  least  square  sense.     Portions  of  the 
data  output  for  various  tinnes  are  shown  below  in  Table  I. 
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Table  I.     Portions  of  data  output. 


Time 


{S) 


dp 


192.000 


200.000 


218. 000 


240.000 


15618584 
93124562 
81713672 

14986779 
93515921 
81496110 

. 13411623 
.94365799 
. 80819492 

. 11553407 
.95309533 
. 79952632 


1.2487292    0.00000196 


1.2493923    0.00000038 


0.0008078 


1.2494577    0.00000001    0.0000527 


1.2496626    0.00000003    0.0001062 


0.0003095 


Time  is  in  seconds,    and  lengths  are  in  earth  radii. 


The  smallest  error  sphere  corresponds  to  time  200.  000.       The 

radius  of  this  sphere  is  found  using  (6).         Here 

1 

/    0.00000001    \^  „     nnnnc->-7  *v  J- 

dP    =    ( ^ =•  )         =  0.  0000527    earth  radii, 

which  is  about  0.  2  miles. 

Since  time  was  used  as  the  parameter,    the  total  computed  change 
in  parameters  at  each  station  is  an  indication  of  the  timing  error  between 
stations.     The  variations  in  time  may  be  seen  below  in  Table  II. 


Table  II.     Adjusted  parameters  (time). 


Time 


Manhattan 


Menlo  Park 


Oakley 


St.    Louis 


192.000 
200. 000 

208. 000 
218. 000 

230. 000 
234. 000 

240.000 
242.000 


192.  000 
200. 000 

208. 000 
218. 000 

230. 000 
234. 000 

240. 000 
242.  000 


192.  633 
199.556 

207.965 
217.816 

229. 713 
233.967 

239.845 
242.536 


192.  514 
199.653 

208. 140 
218.073 

229.994 
234.085 

239.481 
241.476 


235.055 

239.928 
240.864 


EFFECTS  OF   VARIATION  IN  MEASURED  ORBIT   PARAMETERS 


'  Using  the  previously  computed  sequence  of  satellite  position  vectors 

and  corresponding  error  spheres,    the  effect  of  variations  of  the  measured 
parameters,      P     and    9  ,    on  sonne  of  the  other  orbit  parameters  nnay  be 
analyzed. 

The  sequence  of  satellite  vectors  is  first  corrected  for  the  rota- 
tion of  the  earth  during  the  elapsing  times  of  observation,    which  will 


affect  an  instantaneous  set  of  position  vectors.     The  vectors  are  cor- 
rected  using  the  transformation 


fi 


rcos  U). 


sin  (I). 


■sin  OJ. 
1 

cos  a>. 

1 


(7) 


where    (0.    =    n  x  (t.  -  tj  )     and    fi     is  the  angular  velocity  of  the  earth 
in  radians  per  second.      (_x,    y,    z  J^    is  the  satellite  position  vector 

at  time    t. . 

1 

The  determination  of  the  best  orbit  plane  was  calculated 
using  the  least  square  technique  described  by  Kirmser  and 
Wakabayashi  ^zl.       The  spherical  co-ordinates  of  a  unit  normal  to 
the  orbit  plane  are 

e  '    =     sin'^    n 


m 


and 


where 


4.1    =    tan>    l-^ 


n 


and 


U-r^v  -  (e-e)(  Ti-  n) 


n 


n/    i  +  (^)^+(4:)' 


'n  '        *  m' 

The  position  vectors  determined  in  (7)  are  transformed  into 
the  orbit  plane  co-ordinate  system  with  the  transformation 


0  0      \    /-sin  <j)  •      cos  <j)  '    0\      f  i 

sin  e  '        cos  e  '1/  -cos  (|)  '    -sin  (|)  '     oj    j  ^ 

■  cos    0'        sin  e  '/ \       0  0  l/     I    C 


The  basis  of  this  new  co-ordinate  system  is  the  ascending 
node  vector,  the  normal  to  the  orbit  plane  (orbit  vector),  and  their 
mutual  perpendicular. 

Geometrical  Approach 

The  effect  of  variations  in  the  measurable  parameters  may  be 
considered  through  geometry  alone.     The  equation  of  a  focal  ellipse  in 
polar  co-ordinates,    Fig.    3,    is 


where 


P    =    e  (p  +  P  cos  (  e  -  co)  )  (8) 


P  =  radius  vector, 

e  =  eccentricity,  •     - 

9  =  angle  with  respect  to  the   ^*   axis, 

p  =  distance  from  the  focus  to  the  directrix, 

and  CO  =  angle  from  the  ^*   aixis  to  the  major  axis. 
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Fig.    3.     Focal  ellipse  in  the  orbit  plane. 

Equation  (8)  may  be  rewritten  as 

ep  +  e  P  cos  (e    -  co)   -  P  =  0  =  f(e,    p,    P,  6  ,    uj). 
Taking  the  derivative  of  this  expression  gives  the  matrix  equation 

\   df_    di     ,      (^  ^l  ,      di  df  d{  ^^^ 

Lap    3rJ    (. 


de 


'    i-^     -  JP    -  K7j      i^P 


(9) 


where 


^^     =     e  cos  (e  -  w  )  -   1,      44-    =  -  e  P  sin  (e    -  w  ), 


•g—    =    p  +  P  cos  ^  e  -  CO  /,     T—    =  e 


de 
f_ 


11 


and 

di 


b> 


=  e  P  sin  (  9  -  o)  ) 


Since  it  is  necessary  to  invert  the  coefficient  matrix  on  the  right 
hand  side  of  (9)  to  solve  for  the  unknown  variations,  three  distinct  sets 
of  data  are  used  to  form  a  square  coefficient  matrix.     Then  (9)  becornes 


(b,-Ma,     f^el     =       I2 


(9-) 


do) 
where    B        is  a  3  x  3  matrix,    and    A    is  a  3  x  2  matrix. 

Since  the  experimental  data  covered  only  a  small  portion  of  the 
orbit  the  columns  of    B    are  almost  the  same,    element  by  elennent, 
and    B    tends  to  be  singular.      Therefore  the  three  distinct  sets  of 
data  used  in  formiing  the  matrix  equation  (9  '  )  are  picked  from  the 
first,    middle,    and  last  sets  of  sighting  data  to  reduce  this  tendency 
toward  singularity. 

However,    the  matrix    B    is  still  very  nearly  singular,    and  for 
small  changes  in  P  and  9,   (9 '  )  gives  quite  large  values  of    de,    dp, 
and  do) .       Because  of  these  large  scale  changes  in  the  orbit  para- 
meters due  to  small  measurement  errors,    the  orbit  parameters  can- 
not be  satisfactorily  deternnined  by  geometry  alone. 
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Dynamical  Approach 


Dynamical  principles  may  also  be  used  in  considering  the  para- 
meter variations.     Parabolas  are  passed  through  three  well- spaced 
points    (I*   ,   t\*  ).     The  parabolas  are 


r.    ♦      - 


a  I    +  a  z    T.    +    03    T^ 

^1     +     02     T.     +     133     T.2 


(10) 


where        T.    =    t.  -  ti 
1  1        ' 


The  radius,    velocity,    and  angles  shown  in  Fig.    4  are  given 
by  the  set 

I    P,  I    =  -I  c^i'    +    ft'  .       . 


and 


I    Pi  I     =  ^/  az'    +    ^i  . 


tan-1     § 
«i 


1        e  sin(cL)-  6,) 


(11) 
(12) 
(13) 
(14) 


Fig.    4.     Vectors  in  the  orbit  plane. 
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Since  in  central  force  motion  the  angular  momentum  is  a 
constant,   we  have 

C    =     I    P,    XP,   [     =     I    P,  (    |p,  I  cos  <}.,.  (15) 

Also 

2E    =     |P,   |2    .-^  (16) 

where    E    is  the  total  energy  of  the  system  and    K    is  the 
gravitational  constant. 

With  these  expressions  some  of  the  orbit  parameters 
are  given  by 

(17) 

-'     -P'Z  (18) 

From  (15)  and  ( 1  7) 

d  (ep)    =    -^—    dC 

and  dC    =     |  P,|   cos  (j),  dP,  +   |  Pi  [   cos  <<>,  dP,  -   |  pj  [  P,  |  sin  <t),  d(}),. 

Now  as        Pi     =    oti     +    j3/    , 

I  Pi  I   dp,    =    a2  da2+  ^2  d/Sz- 

Also  ^.      =      Q-i+Q-jT.   +   ^3    7.^, 

1  '  1  "^      1      ' 

SO  that       4  .   i  Tj      -    o-i 


• 

C2 

^P    =     K 

and 

e    =  ^y  1  +    2^C 

k2 

2 


2dP/ 


>2- 


and  de  ,  =    da  2    =    — ^^^^^ —     =  d^ 

This  gives  -,/        .    o    \    ,„ 

dp      -  2(Q'2  +  g2)  dP, 

tlri      —  ——, — 5 — I — > — r 

(t2-   tj     |P,[ 

From  (14) 

"^^^    "      (T°-  e'cos  (o)  -  9))'     {   8in(a>  -e,)  de+(e  cos  (o)  -  e,).  e2)da) 

+  (-e  cos  (  w  -  e,)  +  e^  de,}  . 


Combining  the  above  gives  the  matrix  equation 
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where 


\   'n  X     -L  2  Pi   cos  <^i  {az  +  /Bz) 

P,  cos  <{>  1  +   ; — ; ;; 

•-  P,(tE-  ti  ) 


m  (e  cos  (9  I  -  0) )   -  eM 


/a 


m    = 


Pi  P I     sin  <j)  1  cos^  <))  1 
(1  -  e  cos  (9  J  -  w)  )2    . 


(19) 


Once  again  it  is  necessary  to  invert  the  coefficient  matrix 
of  the  unknown  orbit  parameter  variations.     Since  it  is  a  1  x  3 
matrix,    three  sets  of  data  must  be  considered  to  form  a  3  x  3 
matrix,    or  a  sqxiare  matrix,    which  may  be  inverted.     (As  before, 
well- spaced  sets  of  data  are  used  since  the  coefficient  matrix 
tends  to  be  singular.)     This  gives  the  matrix  equation 


(A) 


from  which 


=    (B) 


/::) 


=  (A)-1    (B)/l7 


(19') 


(20) 
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The  coefficient  matrix    A    is  nearly  singular,    as  it  was  in 
the  geometrical  investigation.     Although  not  as  bad  as  in  the  first 
investigation,    it  nevertheless  showed  that  for  small  variations  in 
the  measured  parameters     P     and     9  ,     there  were  large  changes 
in  the  other  orbit  parameters. 

LIMITATIONS   ON   ERROR   SPHERE   SIZE 

From  the  previous  work,    the  size  of  the  maximum  allowable 
error  sphere  for  specific  error  tolerances  in  the  various  orbit 
parameters  may  be  calculated. 

From  (l9)  it  is  seen  that  the  coefficient  matrix  of  the  measur- 
able parameters  is  a  1  x  2  matrix.      Two  sets  of  distinct  data  are 
considered  so  that  this  matrix  is  built-up  into  a  2  x  Z  matrix.     Then 
it  may  be  inverted  and  (l9)  is  written  as 


:dP 

de 


•d{ep) 
=    (B)-l    (A)      Jde        /  (21) 


do) 

For  an  error  allowance  of  0.  1  per  cent  in  the  orbit  parameters 
ep,     e,    and    to     it  is  found  from  (2l)  that  the  error  sphere  will  have 
a  radius  of  approximately  1  x  lO"      earth  radii.     The  smallest  error 
sphere  computed  from  the  set  of  experimental  data  had  a  radius  of 
5  X  10      earth  radii  or  0.  2  miles,    which  is  about  five  times  larger 
than  the  error  sphere  required  for  the  specified  error  allowance. 

Typical  values  of  ep,    e,    and    cj     were  obtained  from  sample 
calculations  in  l.2j. 
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The  necessary  increase  in  accuracy  could  perhaps  be  attained 
from  an  analysis  of  the  techniques  used  in  obtaining  and  reducing 
the  experimental  data  as  detailed  by  Kirmser,    Wakabayashi,    and 
Creech  I_3  J  . 

.  DISCUSSION 

This  report  shows  that  the  given  sighting  data  did  not  include 
enough  of  the  orbit  to  determine  the  orbit  parameters,    but  that  the 
various  satellite  positions  could  be  located  quite  accurately.     If 
data  were  available  from  another  set  of  observation  stations,    located 
at  least  one  third  of  the  way  around  the  earth  from  central  USA,    the 
orbit  parameters  could  be  satisfactorily  determined.     This  second  set 
of  stations  would  remove  the  singularity  of  the  coefficient  matrix  that 
"blew  up"  both  the  geometrical  and  dynamical  analyses.     It  would  be 
best  to  have  three  sets  of  stations  located  at  the  third  points  of  the  orbit. 
The  sighting  data  used  in  this  report  covered  less  than  5     (central  angle) 
of  the  satellite  orbit. 

Several  other  points  should  be  made.     The  eccentricity  of  this 
particular  orbit  was  in  the  neighborhood  of  0.  0  2  which  is  almost  a 
circle.     For  a  circle    p  =  oo ,    hence  with  an  orbit  approximating  a  circle, 
p    will  be  quite  large.     Also,    as  an  orbit  approaches  a  circle  the  deter- 
mination of    o)     becomes  meaningless,    for  any  cixis  of  a  circle  is  a 
major  axis. 
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As  stated  earlier,    time  was  used  as  the  parameter  to  determine 
the  timing  error  between  stations.     However,    the  timing  error  values 
must  be  viewed  with  caution,    for  this  timing  error  will  include  other 
errors  such  as  errors  in  the  location  of  the  observer's  positions, 
observational  errors,    computational  errors  because  of  the  parameter 
technique  used,    and  errors  in  the  reduction  of  photographic  glass  plate 
negatives     [_  3  J  .  ,  ' 
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APPENDICES 


2fi 


APPENDIX  I 

The  appendix  contains  the  IBM  1620  programs,    written  in 
FORTRAN  II  language,    and  the  corresponding  flow  charts,   which 
were  used  to  convert  the  experimental  data  to  satellite  position  vectors 
as  explained  on  pp  1-7  of  this  report. 
The  progranns  are  - 

I,       Program  for  normalization  of  the  experi- 
mental sighting  vectors  and  calculation  of  the 
elements  of  the  3x3  matrix,    defined  in  (l'), 
for  each  normalized  sighting  vecotr.     The  in- 
put is  the  set  of  sighting  vectors,    their  cor- 
responding times  of  observation,    and  the  station 
^  from  which  they  were  observed.      The  stations,    or 

observer  locations,    were  numbered  1  through  N 
with  the  number  N  being  given  to  the  reference 
station.     The  output  are  the  elements  of  the  M 
matrix.     This  output  makes  up  a  portion  of  the 
input  for  program  II. 
II.         Program  for  solving  the  system  of  nonlinear  equa- 
tions (4)  by  an  iterative  technique.      The  input  in- 
cludes the  number  of  stations,    the  initial  estimate 
of  the  parameter  (time),    the  observer's  position 
vectors,    and  the  output  of  program  I.     A  typical  set 
of  input  data  is  shown  after  the  listing  of  program  II. 
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Since  derivatives  of  the  M  matrix  are  required,    a  parabola  is 
passed  through  the  corresponding  elements  of  three  M  matrices  to 
form  a  new  matrix  whose  elements  are  functions  of  time.      There- 
fore it  is  necessary  to  read  in  three  M  matrices  for  each  partici- 
pating station.     The  three  matrices  chosen  should  bracket  the  initial 
parameter  (time)   estimate.       For  the  M  matrix  input  the  cards  must 
be  in  the  following  specified  order: 

1.  The  row  one  cards,  from  each  of  the 
three  M  matrices,  arranged  from  lowest 
to  highest  time. 

2.  The  row  two  cards  from  each  of  the 
three  M  matrices. 

3.  The  row  three  cards  from  each  of 
the  three  M  nnatrices. 

This  is  done  for  station  one,    then  two,    and  up  to  station  number  N 
so  that  we  have  nine  cards  for  station  one,    then  nine  cards  for  station 
two,    etc. 

The  output  of  program  II  includes  satellite  position  vectors,    the 
final  value  of  the  parameter  for  each  observer,   and  the  radius  of  the 
error  sphere. 


FLOW   CHARTS 
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PROGRAM   I 
[  Start\ 


Input 


Computations 


Output 


PROGRAM   II 
(Start  J 


Input 


Computations 


h 


K 


1x10 


■8\  >   0 


K=  1,  2,  .  .  .,/N-l 


<   0 


Output 


(End  J 
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C  N  0  r<  i'l A  L  i  Z  A  T  1  C  fvl       D h  1 1 K  I-i  I  N A  T  I  C  N    C  F    ri  A_  1_  k_I  C_E S  _ 

1    Di;/EMSICN    A:-'i(3»3) 

4"Tcr.\'/^Ar'rirr4V5V""Fsv3. 12  > 

Read  in  Sighting  Vectors  (X,  Y,  Z),  .Time_CT)»_and^ 

Corresponding   Station  (ISTA) 


"5  KtAD'  zr»x^  Y'»z  n  n  s  r  A 

Normalization  of  the  Sighting  Vectors 


SCJRTF     (X*X  +  Y*Y  +  Z^<^Z) 


D=     X/A 
C  =  Y/A 

'D'=Z"/A"" 


-Computation  of  the -3x3^ M  Matrix 

AKI  1,1)  =l,u-b^^b 


An( 1 »2 )--o*C 

A  !•  I  (    1   ,  3  )  =^  —  O  ><•  LV 

"MT'r(  z»Tr="Ar-rriT2T' 

_AiM(^»2.)^l.o-C-"C 
"AT^i  (2  »3")  =-C-*i-L' 
AM  (  J»  1  )  =Ai'i(  1,3) 


Ai-;(  3,2  )=AM(2  ,3) 

_^!liA!ALf  I'ii7i2^>i 

13  c   9      l  =  U    3  7" 

-Read  out  the  M-  Matrix  (l — -Row  Number) 


9    P  UKCri    1  ^  ,    AV,  {  I  ,  1  )  »    A  ,v,  (  I  ,  2  )  .    AM  (  I  ♦  3  )  »    T  ,     I  51  A  ♦      _  I 


I  (J    F  CR  i'l A  T     (  3  r  1 4  .  1  G  »  r  6  .  3  >  2  I  2     ) 
Return  to- 5-iar-a- New-Set  of  Sighting  Vectors 


GC    TO    3 . 

' 

3o    tND 
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PROGRAM  II 


*16    5    MAIN    PR'JGUf^!^. 

i^]i^'tNS_I_C;N___j^AJ^  (_j.)_»___CMA  (  3  »  3  I  4  )  >_  CiMo  (  ^  »  ^  »  4  )  »     Cm^  (  3  »  3  *±) 

Ol  MiiNb  i  l;N    T  (  b  )  »    AM  (  3  »  3  »  4  ) '/  "AMPT3V3V4  )  ,   ~A|V,pp  (  3  .  j  » 4  )  ' 
DIMhNSlCN    AI(3»3)»     S55(3»l»4)»    oS(3»i)»    b(:>»l) 


DlMcNSICN     bDb(3.1»4)»     dtiDD  (  j  »  1  ^4  )  »     d(4)»     AA(3»3»4) 

i^I  ■'lLNSJ-.?^LPJ_^A^J_»__AiiLl_»  3_»_4J_»__A_oD  (  4  *J*)  ,  __A  t  L4^»  3  )_ 

DIMENSION    P(  3.3»4)  jAOj^S)  »bD(  f»3V4")  ♦A'"AM(r,3»4) 

Ill!lP^i5Jj5'^L_L&lAl_».QT_(.4_l___,E_A13_.J_»_4_l __. 

niME~N<^,ICMx]9)  »C7(9),Y(9)»C(9)»U(9).h{9),F(9) 

?     Frf^MAT     (F5.3) 

-?    FC^MAT     (3F    14.8) 
_.J^-'^?P'^'AJ  _'^F;4._ipj_F8._3,2_I2)_ 

__iA_f_?il'^A  L_Li_4_»__t_o  ?>  1 

26   FClRMAf     {Ye  16.  lb*    312  y       ' 

Read  in  Number-ofStations-and-Initial  Parameter  Estimates 

2_91^'^.^AQJAy^}^JJ5_} 

Read  in  Step  Factor 


AM  =  N 


Read  in  Observer  Position  Vectors  (K-Station  Nximber) 


114    READ    3,    P(i.l»K)»    P(i,2»K.),    P(1.3»(^) 

K  U  =  c 

205    DC    6    1=1,3 

Read  in  M  Matrix  (Program  I  Output) 


.A__RF_vU__4j>._X_LJJ.i_J(J_I+_l)  ,  _  X  (  L+6_L!lT  LLJ.»_.15-I^.  -_L, 

Ml.)     ^NU    +     1 
DC     13     1=1  .7.3 


Pass  Parabola  Through  Corresponding  Elements  of  M  Matrix, 
Then  M  Becomes  a~ Function  of  Time 


._Ui-IJ_=XJ_D^J_LJ_LI_Lr_LlJiJ.L1^XT_(_IJ_rJ_C3_LJJ 

H(l)=X(i+l)/((T(2)-T(l))*(T(2)-T(3))) 
K  1  )=X(  i+2  )/(  (T(31-T(  1  1  HW  T(  Hl-TOi  1  1 


G(  I  )=T(2  )i-T(  3)*U(  I  )+  f  (  1  )*T(  3)*n(  I  )+T  (1  )irT(2  )*h  (  I  ) 

-tJ_LJ_=r_(.i_IiZl_+_LJ_'iJ_Lr_UJJ.J_+iilJLi_+_ll^J_l2^JllJ.l_+_lJ_iXli-_Ii^j  i^i  i_IJ  I 

1 3    C(  I  )=U( I )+n( I )+F{ 1 ) 

iir._J^_±EJLi2. 


-26. 


MU    =    3^^1-2 

CMA     (L>  I  »  ISTA)     =    GdMu) 


CMd     (L»i»ISTA)  =     y(ML) 

1_2    CMC_J_L_»Jj.l_5TAl  =__CiMU_)_ 

rF"(MU-3*T!T  "  20^,213,215 
2]  5    DC    Jfi    K  =  1,M 


Set  All  Parameters  Equal  to  Initial  Estimate 


T^     T  (  K  )  =  T  (  S  ) 

T6"nc'r7"K  =  i  .n"  " 

f)?.    17    l_=l_,3 

DC     17    J=l,^ 

Expressions^for-M^   M  V-andJML' as-JFunctions  of  Time. 


AM  [i_»_J_'_K_)  =CM  A_[  I_,  J_,  Lll'i-J-'LD  ii_',::'_»_5_l?'J_iAl±^i''i^-'_L»_-:>_»>.  )JlIJ_'^J_^rJ_ 

AMPn  ,J^<r=CMb"(  r",J»K  )+2".*(  Ci*iC(  L'*J*K)-kjC^)) 
1_7_  AMPP(  I_,_J_,_K_)_=_2_.*_(_C_f-^C_(_]  »_J_,lS_)_I 

"DC  "18   1  ="l  ,  ? 
DC    ]R    J=l,3 


A  (  I  ,  J  )  =  "  .  C: 
_DC    JP     '<=l,i^I 

Computation  of     A  and    A 


18    A(  I  ,J)=A(  I  ,J)+AIvi(  I  ,J,K) 

AD=(A(l,l)^A(2,^)-'^A(j,3)+A(l,2)*A(^>5)^A(3,l) 


1  +  A(  1,3)*A(  3,2  )*A(<i  ,1  )-A(  1,  j)^A(^i,2  )-A(  3,1  ) -A  (  2  ,  3  ) 
.2_*A  L3_,  2_).*AiJ_5JJr^AJ_3j_3Jj''i:t\J_^_,_iJ.*A(  1,»  2_1J 

AI  ( 1,1 )  =  { A ( 2,2 )*A ( 3,3 )-A( ^ , J )^A( 3,^  )  ) /Au 
.__AJ_(_l_»?J_r-iAiJ.»ZliLAr3_,_'lJ.-AJ_l_,_3_)LlrA13_».^_LlV_H^ 

AI  (  1,3)  =  (  A(  i  ,2  l*'^!  2,3)-A(  1  ,3  )  "-A(2,2)  ) /Au 

A  I  (  2  ,  1.  )  =-  (  A  (  ?  ,  1  1  -^A  (  3  ,  3  )  -A  (  2  ,  3  )  ^-A  (  3  ,  1)  )  /mJ 


A  I  (  2  ,  2  )  =     (^(1,1  )  -"-A  (  3  ,  3  )  -A  (  1  ,  3  )  ^i  A  (  3  , 1  )  )  / AD 

J^J12_,_3j_f-_LA_Li_»Ll*AL?»J^lrAJJ_j_3j_«:M_2_».iJ_I/Ap__. 

AI  (3,1  )=     (  A(2,l  )*A~(  3,2)-A(  2,2  )*A(  3,1)  ) /AD 

A  I  (  3 , 2  )  -r_[ALJ_iL'  *AL3_ »iLir_^Lli_?._)_*^-l'?AD_l/AlI_ 

AI(3,3)=     {A(l,l)     *A(2,2  )-A(  1  ,2  )*A(  2,1  )  )7aiJ 


-CompxiLaLion  uf    (M^;) — {-R-T")- 


OOJ.Z_^J__K^Ul^_ 

UC^^O     i -i ,3 
^SSJ  l_*j_:>K  )_=0_*_G 

2?n    5SS(  I  ,  I  ,K)=SS5(  I  ,1  ,K)+AM(  I  ,  J  ,  K  )  »^- (  1  ,  J  ,  iC  ) 


Computation  of    2(M^  )   {  Pt<.  } 


k^i_-^:js._ 


DC    25    1=1,^ 
_S^LLi:_lrII.J3 

DC    •'^    <  =  T,!\ 
25    SS(  I  ,  1  )^bo(  I  ,  1  )  -r5SS(  1  ,1  ,K,) 
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Computation  of   S  =  A"^  S(Mj^)  {  Pj^  } 


uC    :)0    1-1  »3 
__  ^  (  I  *Xl=^*'l ' 

Computation  of      P„     -  S      I 

DC    "^^    K=l  ,N 
jDC._3i3.^J_=l.».3. 

35    bB(  i  »J»K)  =P(  l.J»i<)-S(  J»l  ) 


Computation  of  "Satellite  Distance 

31     5DS    =    0.0 


DC    32       1=1,3 
J_a_JiPJi..iL  ^  S  OS +5J_  J_»_ -Ll*iit  J_ » 1-L 
SD    =    S0KTF(SD5) 


"Computation  of  t^^  -  ^} 

is. 


no    40    <=! ,M 


DC    40     1=1,3 

_4_c__3B  I  (  u  l.»K_L=.PJA_'_l_'i^l--'?_U  vLL— 
Computation  of  (M^,*    )   {  P^,  -  S} 


K    '    '  *  K 

DC    45'"K=r,NT 
DC    45    1=1,3 


BBBri( I ,1 ,K)=0.0 

_DC._4^_J=.L,3. 

4  5    HBBR(  I  ,1  ,<)=BBRB(  I  ,1  ,K  )+AMP(.I  ,  J  ,K  )  *dBB  ( J  ♦  1  ,K  ) 

rp 

<E;omputation  of-  [  Pj^     -  Sj  {M^^)  {  P^^  *  S  > 


_D.C_„5-0-)CjlL»W-L.. 


P  ( ,<  )  =  0.  0 

j[)C._50_J=r,.3- 

50    B(K)=B(K)+DD(1,I  ,i<)*BBBB(  I  ,  1  ♦  K. ) 


IZIomputation  of  ~2(M;j-*  )    A"    " 


-UC-   5  5    K  =  ]  .N-L 


DC    55    1=1,3 
JiC._5-5._^ri,i 

AA(  I  ,J,K) =0.0 

-DC-_55-i_=JL»3- 

55    AA(  I  ,J,K)=AA(  I  ,J,K)+2.*{AN!lP{  I  ,1    ,  K  )  ^^  A  I  (  L  ,  J  )  ) 


Computation  of  [P   '^  -  S^  J(2Mj^'   A"^j 

J2C__6_0_jLrJLii^I 

DC    60    J=l,3 
Ai^ALUAiAl=Jdj'J^ 

DC    60    1=1,3 
fiO    AAAd  ,J,.<)=AAA(  1,J,K)+BB(  1  ,  I  ,i<)  ^^AA  U  ,  J  ,,-s  1 
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^-^T  ,  ..V,   T-A-1 


Compulation  oill'Pj^'-S"^  J  C2Mj^^X'*)(Mj^j^  )   {  Pj^^  -  S}" 


DC    65    K=1,NI 
DC    6_5_KK=1  ,_N_I_ 

"d ("<"♦<<)  =0.0'  " 

_    DC    65     1=1 ,^ 
"65'  D(tC  ♦<kV=D(  K.  KkV+AAATi  »  I  ,<)*"RRBB(  I  »1»X.K) 

K     -  S    -I   ^"-   K 

^-.Zii_5!l»N.i 

DC"7l)  j=T,y" 

A b  (  1  ♦_J_»_l<.j  =0_«P_ 

PC^7C    1  =  1",  3"  '       " 
70    AB(l»J»K)=AH(l»J>K)+Db(]»I>K  )^A;-;FP(  I  »J>i<) 

Computation  of  (_P   "^    "  ^''' J  t^'  ^    ^  ^KK  "  ^^ 
'dC    75    <  =  1,NI 
DC_75_J<<iLljJ^J 

ABB (K»KK) =0.0 
DC_i&_UJLjL3 


Computation  of-[_P— ^-=-5    J  ( M '^.r )- 


7  5    Abb(K»,<.K.)  =ABb(,';,^K)+Ab  (1,1  ,^)-"-bbH(  i  »i,KK  ) 

OC,__9_0_><=±iM *^- 

DC    90    KK=1»NI 


p  r-    A  P  (  <  ,  1^  <  )  =  n  (  K  ♦  !^  <  ) 

'^^_Tj^_'iv 

PS    AP(i<,<K  )=0(K»K<  )-AaB{K»KK) 

.9o__cy^J.lij'k:; 

Computation  of    E  =  S  f  P^,"^  -  S^^J  (M^  )   {  P.,  -  S} 


95    DC    iOO       f;^i,i\ 
j5C__iO.O__J_=J_»3_ 

"  EA(  I  ,r,  K)     =    u.O  ~ 
0_OAOJJ___J_=1_\3_ 

100    EA(I,1,<)     =    EA(  I  ,1  »<)+AM(  I  ,J,K)^Bbb(  J,l  ,f<) 
DC    105       K=1»M 


•^  P  (  <)     =    0  .  n 
DC  J_Q5l l=X}.l 

105    FB(K)     =    FB( K  )+BR( 1  ,  I  ,K)*EA (  I ,1,K) 
E_=_Jls,^ 

DC     110       ,<-l,M 
-m)_L_^^E  +  t:B(KJ 


Solution  of  N-  1  Equations    (a..){A|}   =   {^} 

DC    135    "1  =  1  tNI 
,13_5__AEAJ_j.N_Lr3i_IJ 

DC    125    ^=l,i\il 
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^0    ni     1  =  1  ^NT 


DC    111    J=l .N 
_1L1_  j':?.!  J_  >  J  _)__=__Af  J_l,  J^L_- 

DC    115    I  =  l,iNI 

_y.(_l.»_NIJ__r-^«_t^-(_L»A-L 

115    D( i  »K)     =    AE(  I ,NI ) 

DC    120    M=2,NI 


DC    120    I=iV,NI 

ilC^JZQ_JfL^:^ii 

120    D(I»J)     -    D(M-1  »;'i-l)*D(  I  »J)     -    D(  i  ».v,-i)-><D(;'i-l  O) 

Change  in  Parameter  for  Each    Station  {a  i)- 


12^    DT (K)     =    D(NI «N)/D(NI >NI ) 


]25    rCiVTINUE 

N  C_„r.^N;<+ 1 

DC    199    K=1»NI 


Comparison  of    Aip.  with  a  small  number:"  If'.all  A^  „    .are  smaller 
than  the  number,    process  output  and  return  to  the  start  of  the  pro- 
gram  for  a  new  set  of  data.     If  not7"c6ntTnue  iteration  using  the 
newly  adjusted  parameters. 

IF     (ARSF(DT(  K)  )-1.0E-08)     19-/,1^9»17& 
_19J5__CCWIINiiE 

EP    =    ARc.F     (E) 


159    PRINT    167,     T( 5) 
.16JCl_-PRJ.^LI_.lbi_,__c 

161    FCRMAT     (iOh    EkkCR    IS    Li4.o). 
_162:__PR-LNl_j._(23_»_Jii-_i_^D ^ 

163  FCRMAT     (  Yri    TKiAL    12  »    ^h    SO    cll^.o) 

164  FCRf/AT     (jEi4.o>     Fci«3) 


PUNCH    165,     T ( 5) ,    E,    bD,    UR 

?_yNCH_J.64_,_.SiJ„»i_ljL_3.r2_rI.l_,__S_L3,,_l±,_J_L3J_ 

PUNCH    167,     (DT(i<)  ,     K  =  l  ,NI) 

nC__16.8.J^5__J-,L!^! 

16R    PUNCH    16°     ,     T(K ) ,     K 

16^     FCRVAT     (4F14.n) 


167    FCRN'AT     (tl4.d) 

_l6-9__F.CJi^lAT__ii^lit..o_i__I^L 

GC    TC    201 
-17_Q-J^C____L7ii.._Jl_=_i_5ui^J 

May  Print  Out  Latest  Set  of  Parameter  Changes 


IF  (SENSE  SWITCH  1  )  171,  175 
.17_I_^ELLIvLI_:L6_Ii__DiIkLi 

175  T(^')=T(K)  +  DEL  ^    DT(K) 

LE__JilC_-^jL__L2_L*__li£ij__i_I6 

175  r-p-L  =T  .C 


3o: 


May  Change  Del,   the  Step  Factor,   to  Improve  Convergence 

177  IF  (.'^EN'^E  SWITCH  2)  178»  200 


178  ACCEPT  2»  OtL 
20_0_GC  T_C__i_6 

end' 


-31. 


TYPICAL  INPUT   FOR  PROGRAM  II 


-7— 2zrOvOtro- 
._^ 


-.4230P268 
-;Tr8F97T4?- 


-.672b730b 
--V7Tr076"5n" 


+  .62  OTTO'^TT 
+.60483236 

■+";6"28'5rsr9" 


.  P  ^  6  fl  7  '>■  7  5  ?  7 
-,  p  6  0  =i  T^  R  9'f>  9^ 
■.25S71  ]  A6f>8 
^.2  5'Tft?A67T' 

•.24523033i:'8 

.?697963?34_ 
-.  24469  5626  6" 

.5544217278 
'.5  5210"84'^99T" 

.5497642085 
%'4''"T4'=''??40  0' 

.^OftfiO] 5304 

•.2846417882 
..■2  8  342^4  5"8l  8" 
-.28?  1753527 
■.98690~5'2"1"'!^9' 
_.9_8  91_52  897  3_ 
.99C46945  7d 
•.082396C01C 

';o75?aT2Tr8"" 

■.071  1Q3073_3_ 

".0'78  3"3T5  304" 

.070830179^ 


■V2^57 
-.2  551 
■  .7  4T2 

.6(^8? 
V6Tr5"8 

.5687 
V4r5B' 
_._4167 

.4302 

.4074 
740786 

.6274 
.6272 
.6"?  70 
.2602 
.258  5' 
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APPENDIX  II 

This  appendix  contains  the  observers'  position  vectors 
and  the  sequence  of  sighting  vectors  and  corresponding  times 
of  observation  from  the  observing  stations.      The  coordinate 
system  has  the    x    and    y    axis  on  the  eqiaatorial  plane  and  the 
z  axis  pointing  north.      The  coordinate  system  is  rigidly  attached 
to  the  earth. 

The  format  of  the  sighting  vectors  is  the  output  of  the 
IBM  650. 
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In  this  report  a  sequence  of  discrete  unsynchronized  satellite 
sighting  vectors  from    N  observation  stations  located  on  the  earth 
was  reduced  to  a  set  of  satellite  position  vectors.      The  effect  of  errors 
in  the  observation  parameters  on  orbit  parameters  is  discussed. 

The  theory  used  leads  to  a  set  of  N- 1  non-linear  algebraic 
equations  in  terms  of  a  set  of  parameters.     This  system  is  solved 
using  a  multi- dimensional  Newton's  method.       Since  the  method  is  of 
an  iterative  nature  a  program  for  computation  was  written  for  the 
IBM   1620. 

The  smallest  error  sphere  obtained  about  a  satellite  position 
vector  was  approximately  0.  2  miles.     It  was  found  that  for  a  variance 
of  0.  1  per  cent  in  orbit  parameters  the  error  sphere  must  be  about 
0.  04  miles,    or  one  fifth  as  large  as  the  smallest  error  sphere  obtained 
from  the  experimental  data. 

Experimental  position  data  obtained  from  segments  of  orbits  sub- 
tending central  angles  of  5° are  ill-conditioned  for  the  computation  of 
orbit  parameters.     The  need  for  super  accurate  data  could  be  avoided 
by  locating  an  additional  pair  of  observation  stations  a  sizable  fraction 
of  the  earth's  circumference  away  from  the  first  pair. 


